Compiled under R version 4.1.2 (2021-11-01)

WARNING: edit the working directory to your preferred folder.

This document details all analyses performed in R for the study:
Legendre, L. J., S. Choi, and J. A. Clarke. The diverse terminology and complex evolution of reptile eggshell microstructure.

For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at .

Loading packages

library(ape)
library(castor)
library(evobiR)
library(ggtree)
library(phytools)
library(RColorBrewer)

Data and tree

tree<-read.nexus("treewhole_newversion.trees.nex")
plotTree(tree, fsize=0.5,lwd=1,type="fan",ftype="i")

data<-read.table("Datawhole_newproject.txt", header=T)
setdiff(tree$tip.label, data$Taxon) # taxa and data match
## character(0)
rownames(data)<-data$Taxon
datanew<-ReorderData(tree, data)

Ancestral state reconstruction: discrete trait (soft/semi-rigid/hard)

thisstudy<-datanew$Type_2021; names(thisstudy)<-rownames(datanew)
norell<-datanew$Type_Norell; names(norell)<-rownames(datanew)
legendre<-datanew$Type_Legendre; names(legendre)<-rownames(datanew)

# Color palette
cols<-setNames(c("royalblue","green3","red3"),c("Soft","Semi-rigid","Hard"))

# Visualize tree with node numbers
ggtree(tree) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

Using character coding defined in this study

  • Perform SIMMAP using phytools – 1000 iterations, using AIC to select best model (code modified from Liam Revell, see here)
x<-thisstudy

aic<-function(logL,k) 2*k-2*logL
aic.w<-function(aic){
    d.aic<-aic-min(aic)
    exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
}

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0010672548  0.0005336274  0.0005336274
## Semi-rigid  0.0005336274 -0.0010672548  0.0005336274
## Soft        0.0005336274  0.0005336274 -0.0010672548
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0010457689  0.0003790695  0.0006666995
## Semi-rigid  0.0003790695 -0.0008898238  0.0005107544
## Soft        0.0006666995  0.0005107544 -0.0011774539
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard  Semi-rigid          Soft
## Hard       0.0000000000  0.00000000  0.0000000000
## Semi-rigid 0.0058104610 -0.00962594  0.0038154791
## Soft       0.0004728781  0.00000000 -0.0004728781
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
logL
##        ER       SYM       ARD 
## -56.06148 -55.75401 -43.99790
AIC<-mapply(aic,logL,c(1,3,6))
AIC
##       ER      SYM      ARD 
## 114.1230 117.5080  99.9958
AIC.W<-aic.w(AIC)
AIC.W
##           ER          SYM          ARD 
## 0.0008548391 0.0001573372 0.9989878237
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
##   1   0 999
o1<-make.simmap(tree,x,model="ARD",nsim=999)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard  Semi-rigid          Soft
## Hard       0.0000000000  0.00000000  0.0000000000
## Semi-rigid 0.0058104610 -0.00962594  0.0038154791
## Soft       0.0004728781  0.00000000 -0.0004728781
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
o2<-make.simmap(tree,x,model="ER",nsim=1)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0010672548  0.0005336274  0.0005336274
## Semi-rigid  0.0005336274 -0.0010672548  0.0005336274
## Soft        0.0005336274  0.0005336274 -0.0010672548
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treesstudy<-c(o1,o2)

objstudy<-describe.simmap(treesstudy)
objstudy
## 1000 trees with a mapped discrete character with states:
##  Hard, Semi-rigid, Soft 
## 
## trees have 18.211 changes between states on average
## 
## changes are of the following types:
##      Hard,Semi-rigid Hard,Soft Semi-rigid,Hard Semi-rigid,Soft Soft,Hard
## x->y           0.003     0.001          10.112           6.631     1.462
##      Soft,Semi-rigid
## x->y           0.002
## 
## mean total time spent in each state is:
##              Hard   Semi-rigid         Soft    total
## raw  8019.4557598 1739.0045051 3095.1775697 12853.64
## prop    0.6239055    0.1352928    0.2408017     1.00
plot(objstudy,type="fan",fsize=0.01,lwd=1,ftype="i", colors=cols,ylim=c(-2,Ntip(tree)),offset=20, part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Here, a semi-rigid eggshell is the ancestral state for almost all major reptile clades – reptiles, lepidosaurs, squamates, turtles, archelosaurs, archosaurs, dinosaurs,ornithischians, and saurischians. However, there are only 6 taxa with semi-rigid eggshells out of 203!

=> Some of these taxa may have a excessive influence on this result – probably the two sauropodomorphs with that eggshell type, since they are the closer in age to the root and to the aforementioned subclades.

Using character coding from Norell et al. (2020)

x<-norell

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013811074  0.0006905537  0.0006905537
## Semi-rigid  0.0006905537 -0.0013811074  0.0006905537
## Soft        0.0006905537  0.0006905537 -0.0013811074
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013651414  0.0004938598  0.0008712815
## Semi-rigid  0.0004938598 -0.0010970016  0.0006031418
## Soft        0.0008712815  0.0006031418 -0.0014744233
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                   Hard   Semi-rigid         Soft
## Hard       0.000000000  0.000000000  0.000000000
## Semi-rigid 0.007423526 -0.012130029  0.004706502
## Soft       0.001380656  0.001222171 -0.002602827
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
logL
##        ER       SYM       ARD 
## -68.56111 -68.13921 -55.04396
AIC<-mapply(aic,logL,c(1,3,6))
AIC
##       ER      SYM      ARD 
## 139.1222 142.2784 122.0879
AIC.W<-aic.w(AIC)
AIC.W
##           ER          SYM          ARD 
## 1.999592e-04 4.126532e-05 9.997588e-01
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##   ER  SYM  ARD 
##    0    0 1000
treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                   Hard   Semi-rigid         Soft
## Hard       0.000000000  0.000000000  0.000000000
## Semi-rigid 0.007423526 -0.012130029  0.004706502
## Soft       0.001380656  0.001222171 -0.002602827
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
objnorell<-describe.simmap(treesnorell)
objnorell
## 1000 trees with a mapped discrete character with states:
##  Hard, Semi-rigid, Soft 
## 
## trees have 23.847 changes between states on average
## 
## changes are of the following types:
##      Hard,Semi-rigid Hard,Soft Semi-rigid,Hard Semi-rigid,Soft Soft,Hard
## x->y               0         0           8.254           5.153     5.562
##      Soft,Semi-rigid
## x->y           4.878
## 
## mean total time spent in each state is:
##              Hard   Semi-rigid         Soft    total
## raw  7730.2217576 1.110435e+03 4012.9805883 12853.64
## prop    0.6014034 8.639076e-02    0.3122058     1.00
plot(objnorell,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Completely different result – let us check which taxa are different from the previous ASR…

datanew[c(which(datanew$Type_2021!=datanew$Type_Norell)),]
##                                     Clade                 Taxon
## Massospondylus        Non-avian_dinosaurs        Massospondylus
## Lufengosaurus         Non-avian_dinosaurs         Lufengosaurus
## Chelonoidis_niger              Chelonians     Chelonoidis_niger
## Testudo_marginata              Chelonians     Testudo_marginata
## Testudo_graeca                 Chelonians        Testudo_graeca
## Stigmochelys_pardalis          Chelonians Stigmochelys_pardalis
## Caretta_caretta                Chelonians       Caretta_caretta
##                       Eggshell_thickness  Egg_mass  Type_2021 Type_Norell
## Massospondylus                       100 158.00000 Semi-rigid        Soft
## Lufengosaurus                         85 500.00000 Semi-rigid        Soft
## Chelonoidis_niger                    400 101.64000       Hard        Soft
## Testudo_marginata                     90  16.63200       Hard        Soft
## Testudo_graeca                       165  16.51650       Hard        Soft
## Stigmochelys_pardalis                210   6.37875       Hard  Semi-rigid
## Caretta_caretta                       80  37.20087       Soft  Semi-rigid
##                       Type_Legendre
## Massospondylus                 Hard
## Lufengosaurus                  Hard
## Chelonoidis_niger              Hard
## Testudo_marginata              Hard
## Testudo_graeca                 Hard
## Stigmochelys_pardalis          Hard
## Caretta_caretta                Soft
# Only 7 taxa are different – two non-avian dinosaurs and five turtles

# Changing character state for the two non-avian dinosaurs
x[c(21,22)]<-"Semi-rigid"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                   Hard  Semi-rigid        Soft
## Hard       -0.00144622  0.00072311  0.00072311
## Semi-rigid  0.00072311 -0.00144622  0.00072311
## Soft        0.00072311  0.00072311 -0.00144622
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013771501  0.0005401113  0.0008370389
## Semi-rigid  0.0005401113 -0.0013929825  0.0008528712
## Soft        0.0008370389  0.0008528712 -0.0016899101
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid         Soft
## Hard       0.0000000000  0.0000000000  0.000000000
## Semi-rigid 0.0068533111 -0.0122551684  0.005401857
## Soft       0.0006039133  0.0006682603 -0.001272174
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##   ER  SYM  ARD 
##    0    0 1000
treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid         Soft
## Hard       0.0000000000  0.0000000000  0.000000000
## Semi-rigid 0.0068533111 -0.0122551684  0.005401857
## Soft       0.0006039133  0.0006682603 -0.001272174
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
objnorell<-describe.simmap(treesnorell)
plot(objnorell,fsize=0.5,lwd=1,ftype="i",colors=cols,ylim=c(-2,Ntip(tree)),part=0.95)
add.simmap.legend(colors=cols,prompt=FALSE,x=20,y=180)

All major clades are now recovered as semi-rigid again => strong influence of the two sauropodomorphs.

Using character coding from Legendre et al. (2020)

x<-legendre

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007636699  0.0003818350  0.0003818350
## Semi-rigid  0.0003818350 -0.0007636699  0.0003818350
## Soft        0.0003818350  0.0003818350 -0.0007636699
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0006269773  0.0000000000  0.0006269773
## Semi-rigid  0.0000000000 -0.0005633914  0.0005633914
## Soft        0.0006269773  0.0005633914 -0.0011903687
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001519306  0.0000000000  0.0001519306
## Semi-rigid  0.0019150212 -0.0019150212  0.0000000000
## Soft        0.0019555075  0.0005401052 -0.0024956127
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
logL
##        ER       SYM       ARD 
## -45.50400 -43.12910 -41.26503
AIC<-mapply(aic,logL,c(1,3,6))
AIC
##       ER      SYM      ARD 
## 93.00800 92.25820 94.53006
AIC.W<-aic.w(AIC)
AIC.W
##        ER       SYM       ARD 
## 0.3422281 0.4978884 0.1598835
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
## 342 498 160
l1<-make.simmap(tree,x,model="ER",nsim=342)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007636699  0.0003818350  0.0003818350
## Semi-rigid  0.0003818350 -0.0007636699  0.0003818350
## Soft        0.0003818350  0.0003818350 -0.0007636699
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
l2<-make.simmap(tree,x,model="SYM",nsim=498)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0006269773  0.0000000000  0.0006269773
## Semi-rigid  0.0000000000 -0.0005633914  0.0005633914
## Soft        0.0006269773  0.0005633914 -0.0011903687
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
l3<-make.simmap(tree,x,model="ARD",nsim=160)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001519306  0.0000000000  0.0001519306
## Semi-rigid  0.0019150212 -0.0019150212  0.0000000000
## Soft        0.0019555075  0.0005401052 -0.0024956127
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treeslegendre<-c(l1,l2,l3)

objlegendre<-describe.simmap(treeslegendre)
objlegendre
## 1000 trees with a mapped discrete character with states:
##  Hard, Semi-rigid, Soft 
## 
## trees have 10.214 changes between states on average
## 
## changes are of the following types:
##      Hard,Semi-rigid Hard,Soft Semi-rigid,Hard Semi-rigid,Soft Soft,Hard
## x->y           0.411     4.143            0.14           0.061     3.772
##      Soft,Semi-rigid
## x->y           1.687
## 
## mean total time spent in each state is:
##              Hard   Semi-rigid         Soft    total
## raw  9104.1361000 220.33623291 3529.1655017 12853.64
## prop    0.7082926   0.01714194    0.2745655     1.00
plot(objlegendre,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Reptiles, lepidosaurs, and squamates are recovered as ancestrally soft-shelled, while turtles, archelosaurs, archosaurs, and all dinosaur clades are recovered as hard-shelled.

The pattern seems to follow whichever character state is coded for the two sauropodomorphs Lufengosaurus and Massospondylus.

To test this hypothesis, we must remove all other taxa (n = 7) that are not coded the same in all three studies, and see if the pattern is still present.

  • Remove all taxa with differences in coding, except the two sauropodomorphs
remove<-rownames(datanew[c(which(datanew$Type_Legendre!=datanew$Type_Norell)),][c(3:9),])
x<-thisstudy[!names(thisstudy)%in%remove]
treenew<-drop.tip(tree, setdiff(tree$tip.label, names(x)))

# Visualize tree with node numbers
ggtree(treenew) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

Replicate analysis on the new tree 3 times, with identical coding for all taxa except the two sauropodomorphs

  • Coded as semi-rigid
logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0009105572  0.0004552786  0.0004552786
## Semi-rigid  0.0004552786 -0.0009105572  0.0004552786
## Soft        0.0004552786  0.0004552786 -0.0009105572
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007204487  0.0000000000  0.0007204487
## Semi-rigid  0.0000000000 -0.0007881215  0.0007881215
## Soft        0.0007204487  0.0007881215 -0.0015085703
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0001356366  0.000000000  0.0001356366
## Semi-rigid  0.0033499438 -0.003349944  0.0000000000
## Soft        0.0016634621  0.000864281 -0.0025277432
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
##  75  97 828
s1<-make.simmap(treenew,x,model="ER",nsim=75)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0009105572  0.0004552786  0.0004552786
## Semi-rigid  0.0004552786 -0.0009105572  0.0004552786
## Soft        0.0004552786  0.0004552786 -0.0009105572
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
s2<-make.simmap(treenew,x,model="SYM",nsim=97)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007204487  0.0000000000  0.0007204487
## Semi-rigid  0.0000000000 -0.0007881215  0.0007881215
## Soft        0.0007204487  0.0007881215 -0.0015085703
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
s3<-make.simmap(treenew,x,model="ARD",nsim=828)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0001356366  0.000000000  0.0001356366
## Semi-rigid  0.0033499438 -0.003349944  0.0000000000
## Soft        0.0016634621  0.000864281 -0.0025277432
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treenewsr<-c(s1,s2,s3)

objsr<-describe.simmap(treenewsr)

plot(objsr,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Most inner clades are now recovered as soft-shelled – most likely an artefact due to the position of Mussaurus – closer to the other inner nodes than any other taxon, since it is the oldest specimen in the sample.

  • Coded as soft
x[c(21,22)]<-"Soft"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008576981  0.0004288491  0.0004288491
## Semi-rigid  0.0004288491 -0.0008576981  0.0004288491
## Soft        0.0004288491  0.0004288491 -0.0008576981
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007184639  0.0000000000  0.0007184639
## Semi-rigid  0.0000000000 -0.0005102417  0.0005102417
## Soft        0.0007184639  0.0005102417 -0.0012287055
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001351033  0.0000000000  0.0001351033
## Semi-rigid  0.0024265091 -0.0024265091  0.0000000000
## Soft        0.0017286691  0.0005520113 -0.0022806804
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
##  47 105 848
so1<-make.simmap(treenew,x,model="ER",nsim=46)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008576981  0.0004288491  0.0004288491
## Semi-rigid  0.0004288491 -0.0008576981  0.0004288491
## Soft        0.0004288491  0.0004288491 -0.0008576981
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
so2<-make.simmap(treenew,x,model="SYM",nsim=106)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007184639  0.0000000000  0.0007184639
## Semi-rigid  0.0000000000 -0.0005102417  0.0005102417
## Soft        0.0007184639  0.0005102417 -0.0012287055
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
so3<-make.simmap(treenew,x,model="ARD",nsim=848)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001351033  0.0000000000  0.0001351033
## Semi-rigid  0.0024265091 -0.0024265091  0.0000000000
## Soft        0.0017286691  0.0005520113 -0.0022806804
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treenews<-c(so1, so2, so3)

objs<-describe.simmap(treenews)

plot(objs,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Same result as with semi-rigid coding, unsurprisingly.

  • Coded as hard
x[c(21,22)]<-"Hard"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007868496  0.0003934248  0.0003934248
## Semi-rigid  0.0003934248 -0.0007868496  0.0003934248
## Soft        0.0003934248  0.0003934248 -0.0007868496
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0006471696  0.000000000  0.0006471696
## Semi-rigid  0.0000000000 -0.000567363  0.0005673630
## Soft        0.0006471696  0.000567363 -0.0012145326
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001501237  0.0000000000  0.0001501237
## Semi-rigid  0.0018901522 -0.0018901522  0.0000000000
## Soft        0.0020179175  0.0005399515 -0.0025578690
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
## 354 496 150
h1<-make.simmap(treenew,x,model="ER",nsim=354)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007868496  0.0003934248  0.0003934248
## Semi-rigid  0.0003934248 -0.0007868496  0.0003934248
## Soft        0.0003934248  0.0003934248 -0.0007868496
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
h2<-make.simmap(treenew,x,model="SYM",nsim=496)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0006471696  0.000000000  0.0006471696
## Semi-rigid  0.0000000000 -0.000567363  0.0005673630
## Soft        0.0006471696  0.000567363 -0.0012145326
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
h3<-make.simmap(treenew,x,model="ARD",nsim=150)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001501237  0.0000000000  0.0001501237
## Semi-rigid  0.0018901522 -0.0018901522  0.0000000000
## Soft        0.0020179175  0.0005399515 -0.0025578690
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treenewh<-c(h1, h2, h3)

objh<-describe.simmap(treenewh)

plot(objh,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Ancestral state for dinosaurs and archosaurs, and all major less inclusive clades except pterosaurs, is hard-shelled, as in the Legendre et al. (2020) coding.

To check how strongly are these results influenced by branch length information, we can replicate these analyses using maximum parsimony, which does not consider branch length information, using castor.

  • Coding from this study
studyN<-thisstudy
studyN[studyN=="Soft"]<-1; studyN[studyN=="Semi-rigid"]<-2; studyN[studyN=="Hard"]<-3
studyN<-as.numeric(studyN); names(studyN)<-names(thisstudy)

MPS<-asr_max_parsimony(tree, studyN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPS$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(studyN,sort(unique(studyN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)

MPS2<-MPS$ancestral_likelihoods; rownames(MPS2)<-c(204:369)
  • Coding from Norell et al. (2020)
norellN<-norell
norellN[norellN=="Soft"]<-1; norellN[norellN=="Semi-rigid"]<-2; norellN[norellN=="Hard"]<-3
norellN<-as.numeric(norellN); names(norellN)<-names(norell)

MPN<-asr_max_parsimony(tree, norellN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPN$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(norellN,sort(unique(norellN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=180,fsize=0.8)

MPN2<-MPN$ancestral_likelihoods; rownames(MPN2)<-c(204:369)
  • Coding from Legendre et al. (2020)
legendreN<-legendre
legendreN[legendreN=="Soft"]<-1; legendreN[legendreN=="Semi-rigid"]<-2; legendreN[legendreN=="Hard"]<-3
legendreN<-as.numeric(legendreN); names(legendreN)<-names(legendre)

MPL<-asr_max_parsimony(tree, legendreN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPL$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(legendreN,sort(unique(legendreN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)

MPL2<-MPL$ancestral_likelihoods; rownames(MPL2)<-c(204:369)
  • With the reduced tree, changing only the two sauropodomorphs
xN<-x
xN[xN=="Soft"]<-1; xN[xN=="Semi-rigid"]<-2; xN[xN=="Hard"]<-3
xN<-as.numeric(xN); names(xN)<-names(x)

# Coded as soft
xN[c(21,22)]<-1
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods; rownames(MPx2)<-c(197:355)

# Coded as semi-rigid
xN[c(21,22)]<-2
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods; rownames(MPx2)<-c(197:355)

# Coded as hard
xN[c(21,22)]<-3
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods; rownames(MPx2)<-c(197:355)

Ancestral state reconstruction – continuous character (eggshell thickness)

treedi<-multi2di(tree, random=FALSE)
AET<-log1p(datanew$Eggshell_thickness); names(AET)<-rownames(datanew)
RET<-log1p(datanew$Eggshell_thickness/datanew$Egg_mass); names(RET)<-rownames(datanew)
phylosig(treedi, AET, method="lambda", test=TRUE)
## 
## Phylogenetic signal lambda : 0.999957 
## logL(lambda) : -309.493 
## LR(lambda=0) : 245.84 
## P-value (based on LR test) : 2.09625e-55
phylosig(treedi, RET, method="lambda", test=TRUE)
## 
## Phylogenetic signal lambda : 0.999934 
## logL(lambda) : -228.063 
## LR(lambda=0) : 197.053 
## P-value (based on LR test) : 9.18051e-45

Very strong signal for both absolute and relative eggshell thickness.

phyloplotAET<-contMap(treedi, AET, plot=F)
plot(setMap(phyloplotAET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
     fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
              prompt=FALSE,x=-300,y=280)

# Ancestral states
expm1(fastAnc(treedi, AET))
## Ancestral character estimates using fastAnc:
##         204         205         206         207         208         209 
##   28.079153   28.079153   26.118366   19.854323   36.974371   57.309311 
##         210         211         212         213         214         215 
##   59.388937   38.950551   37.943244   62.261450   63.628751   60.801065 
##         216         217         218         219         220         221 
##   18.564828   17.370830   16.914159   21.363127   14.483029   17.067109 
##         222         223         224         225         226         227 
##   18.167139   27.778883   11.637074    6.394273    5.261709   11.868343 
##         228         229         230         231         232         233 
##   10.874468   12.260022   12.122288   14.011868   16.557668   31.172419 
##         234         235         236         237         238         239 
##   13.895415   14.226503   16.834523   16.374609   14.489393    8.215824 
##         240         241         242         243         244         245 
##    8.889975   16.110137   17.584868   17.641372   11.369380   11.391237 
##         246         247         248         249         250         251 
##   11.091746    7.684418   36.692394  129.109053  132.829774  140.343937 
##         252         253         254         255         256         257 
##  145.689626  141.878765  145.186350   61.134926  159.995521  188.575290 
##         258         259         260         261         262         263 
##  202.108248  218.163318  174.871933  149.316431  250.541379  333.057440 
##         264         265         266         267         268         269 
##  379.854236  544.032218  573.191181  610.483803  607.697486  498.358959 
##         270         271         272         273         274         275 
##   33.924264  333.743365  451.330862  630.020508  397.177198  462.000397 
##         276         277         278         279         280         281 
##  470.274154  478.924007  451.589838   28.662359    2.833103    3.998462 
##         282         283         284         285         286         287 
##   11.270283    3.351885    0.014325   32.566033   69.398144  453.836756 
##         288         289         290         291         292         293 
## 1437.961087 1437.961087  453.836756 1834.612767 1834.612767  453.836756 
##         294         295         296         297         298         299 
## 1523.667987 1544.970306   27.280341   13.898784   72.135564    1.763827 
##         300         301         302         303         304         305 
## 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241 
##         306         307         308         309         310         311 
## 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241 
##         312         313         314         315         316         317 
## 1601.759241 1601.759241  311.981337  473.591727  855.767272  855.767272 
##         318         319         320         321         322         323 
##  855.767272  531.023321  533.631804  533.631804  555.967722 1285.161706 
##         324         325         326         327         328         329 
##  566.272234  605.254074  605.254074  605.254074  605.254074 1101.569284 
##         330         331         332         333         334         335 
## 1391.646422  539.617230  539.617230  539.617230  655.621109  418.158680 
##         336         337         338         339         340         341 
##  418.158680  418.158680  399.605014  399.605014  399.605014  392.933851 
##         342         343         344         345         346         347 
##  399.605014  383.372932  482.275717  482.275717  482.275717  482.275717 
##         348         349         350         351         352         353 
##  482.275717  487.055112  482.275717  672.352457  688.294201  743.214121 
##         354         355         356         357         358         359 
##  749.815848  642.557934 1304.082510  319.361966  778.723541  908.510163 
##         360         361         362         363         364         365 
##  820.749499  423.101767  255.012307  300.635949  343.596583  298.819020 
##         366         367         368         369         370         371 
##  290.372065  222.110312  217.727290  177.850464  166.734252  161.025519 
##         372         373         374         375         376         377 
##  229.672210  234.226968  236.209118  233.541444  232.183846  231.144902 
##         378         379         380         381         382         383 
##  196.457289  216.843344  210.907445  208.757659  210.716135  206.846733 
##         384         385         386         387         388         389 
##  201.747492  251.211818  289.690559  316.985021  253.493514  257.535012 
##         390         391         392         393         394         395 
##  243.362806  151.829900  155.868914  164.019364  144.779678  128.402727 
##         396         397         398         399         400         401 
##  106.958920  100.296111  118.838462  109.666788   73.949794  103.666941 
##         402         403         404         405 
##   99.285261   98.490117  103.797286  101.497666

Ancestral eggshell thickness seems to be intermediate for all major nodes, including archosaurs and dinosaurs. However, many extant lepidosaurs are also recovered as intermediate – likely due to the extremely low values in pterosaurs, reported as lacking a calcareous layer entirely, which shift any thicker eggshell towards the middle of the spectrum.

Furthermore, the pattern strongly follows that of body mass, as expected (e.g. Stein et al., 2019; Legendre and Clarke, 2021).

phyloplotRET<-contMap(treedi, RET, plot=F)
plot(setMap(phyloplotRET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
     fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
              prompt=FALSE,x=-300,y=280)

# Ancestral states
expm1(fastAnc(treedi, RET))
## Ancestral character estimates using fastAnc:
##       204       205       206       207       208       209       210       211 
##  3.409640  3.409640  4.025467  7.755305 20.354753 37.931589 39.819215 25.689681 
##       212       213       214       215       216       217       218       219 
## 39.657045 42.671676 44.850528 53.000367  7.500649  7.357594  7.168380  6.552197 
##       220       221       222       223       224       225       226       227 
##  9.462534  7.377703  1.646676  1.002991  1.610506  0.936549  0.860297  1.828548 
##       228       229       230       231       232       233       234       235 
##  1.596624  2.012481  3.474753  1.682958  2.085547  3.133932  1.180141  1.164824 
##       236       237       238       239       240       241       242       243 
##  7.649844  9.184624 11.645558 12.734058 13.645813 13.020440 16.528768 20.143413 
##       244       245       246       247       248       249       250       251 
## 21.874101 24.369790 23.820456 14.769006  2.599039  6.626096  6.904628  8.138380 
##       252       253       254       255       256       257       258       259 
## 10.904586 10.965487 12.341525  2.047360 11.403766 11.235749 11.089139 11.128995 
##       260       261       262       263       264       265       266       267 
## 10.345142  8.942119 11.414881  8.881152  9.028346 10.830652 10.810594 10.535716 
##       268       269       270       271       272       273       274       275 
## 10.641304  9.256647  2.124890  4.558317  5.340322  6.761045  4.477445  5.802632 
##       276       277       278       279       280       281       282       283 
##  5.664138  5.483004  6.059812  1.692954  0.457460  0.550221  1.153386  0.408094 
##       284       285       286       287       288       289       290       291 
##  0.003995  1.252572  1.164206  1.585486  5.318273  5.318273  1.585486  0.666241 
##       292       293       294       295       296       297       298       299 
##  0.666241  1.585486  1.161423  1.234126  0.894187  0.428783  0.388280  0.130108 
##       300       301       302       303       304       305       306       307 
##  0.766339  0.766339  0.766339  0.766339  0.766339  0.766339  0.766339  0.766339 
##       308       309       310       311       312       313       314       315 
##  0.766339  0.766339  0.766339  0.766339  0.766339  0.766339  2.783888  4.262094 
##       316       317       318       319       320       321       322       323 
##  4.206638  4.206638  4.206638  5.092815  5.314153  5.314153  5.246828  2.804377 
##       324       325       326       327       328       329       330       331 
##  5.263532  4.918171  4.918171  4.918171  4.918171  3.280773  2.041443  5.646772 
##       332       333       334       335       336       337       338       339 
##  5.646772  5.646772  4.638261  8.395933  8.395933  8.395933  8.720703  8.720703 
##       340       341       342       343       344       345       346       347 
##  8.720703  8.804955  8.720703  8.943290  8.062171  8.062171  8.062171  8.062171 
##       348       349       350       351       352       353       354       355 
##  8.062171  8.106581  8.062171  3.124877  3.017738  2.342408  2.259157  2.383106 
##       356       357       358       359       360       361       362       363 
##  0.661660  4.634275  2.112184  1.859717  1.878840  1.176578 10.833970  9.311577 
##       364       365       366       367       368       369       370       371 
##  7.858032  9.853711 10.770694 11.860553 11.582459 13.304214 13.852817 14.336920 
##       372       373       374       375       376       377       378       379 
## 10.818551 10.468502 10.019509  9.677934  9.355059  9.031490 13.464253 11.100662 
##       380       381       382       383       384       385       386       387 
## 10.828725 10.901613 10.572265 10.469195 10.385394  9.574611  8.311332  7.455373 
##       388       389       390       391       392       393       394       395 
##  9.415499  8.811789  9.272832 19.516937 20.275463 19.452284 20.563896 23.256357 
##       396       397       398       399       400       401       402       403 
## 28.370741 30.397420 25.081901 27.482936 41.577907 29.373670 30.524598 32.153403 
##       404       405 
## 31.108659 34.366350

Much lower values for archosaurs and dinosaurs – eggshell thickness seems to have become thinner for a given egg mass at the base of Ornithodira, and later increased in theropods. However, this pattern is dependent on a very small sample of pterosaurs, ornithischians, and sauropods; it is likely that the addition of new specimens attributed to any of these clades would considerably change that pattern.

We can see a strong increase in geckos and in eufalconimorphs – the latter having already been identified in Legendre and Clarke (2021).

References